function s(N)

  close all
  

  %time_steps = [0:N-1]./(N-1)*1.0;


  function hsp = hsourcep(z)
    if z < 0.001
      hsp = 6+(-8).*z+(-7/3).*z.^2+17.*z.^3+(-33/2).*z.^4+(-6).*z.^5+(421/14).*z.^6+(-181/7).*z.^7+(-94/9).*z.^8+(4672/105).*z.^9+(-2507/70).*z.^10+(-139/9).*z.^11+(16360/273).*z.^12+(-21036/455).*z.^13+(-1879/90).*z.^14+(27723/364).*z.^15+(-207577/3640).*z.^16+(-80/3).*z.^17+(183899/1976).*z.^18+(-2356147/34580).*z.^19+(-2293/70).*z.^20;
    else
      hsp = (1/3).*z.^(-1).*(1+z+z.^2).^(-2).*(3.^(1/2).*pi.*(2+z)+3.*z.*(2+z.*(2+z))+(-2).*3.^(1/2).*(2+z).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*(2+z).*log(1+z+z.^2));
    end
  end
  function ksp = ksourcep(z)
    if (z < 0.001)
      ksp = (26/3)+(-7).*z+(-5/3).*z.^2+(28/3).*z.^3+(-107/14).*z.^4+(-7/3).*z.^5+(1327/126).*z.^6+(-571/70).*z.^7+(-26/9).*z.^8+(2413/210).*z.^9+(-7801/910).*z.^10+(-10/3).*z.^11+(5584/455).*z.^12+(-4058/455).*z.^13+(-37/10).*z.^14+(21179/1638).*z.^15+(-637291/69160).*z.^16+(-361/90).*z.^17+(560017/41496).*z.^18+(-7206761/760760).*z.^19+(-899/210).*z.^20;
    else
      ksp = (1/3).*z.^(-3).*(1+z+z.^2).^(-1).*(3.^(1/2).*pi.*(1+z).*((-4)+z.^3)+3.*z.*(8+z.*(8+(-1).*((-6)+z).*z))+(-2).*3.^(1/2).*(1+z).*((-4)+z.^3).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*(1+z).*((-4)+z.^3).*log(1+z+z.^2));
    end
  end


  function resMh = Mh(z)
    resMh = zeros(2,2);
    resMh(1,1) = 1;
    resMh(2,1) = 4;
    resMh(2,2) = z;
  end
  function resh = deqh(z, u)
    h   = u(1);
    hz  = u(2);
    resh = zeros(2,1);
    resh(1) = hz;
    resh(2) = hsourcep(z);
  end

  function resMk = Mk(z)
    resMk = zeros(1,1);
    resMk(1,1) = 2;
  end 
  function resk = deqk(z, u)
    k  = u(1);

    resk = zeros(1,1);
    resk(1) = ksourcep(z);
  end


  z = chebpoints(N, 0, 2);

  
  na2t2  = sa(N, 2);
  %na2t2 = sa_proper(N, 1024, 2);
  na2t4  = sa(N, 4);
  %na2t4 = sa_proper(N, 1024, 4);
  na2t5  = sa(N, 5);
  %na2t5 = sa_proper(N, 1024, 5);

  % ******* h k ******

  reltol = 1e-8;

  uin = [0; 6./4];
  opts = odeset('RelTol', reltol, 'AbsTol', 1e-10, 'MStateDependence', 'none', 'Mass', @Mh);
  [tout, uout] = ode15s(@deqh, z, uin, opts);
  nh2s7 = uout(:, 1);
  nh2s7 = 1/2 + z'.*nh2s7; % with z^2 factored outc  % checked feb27


  uin = [0]; % initial condition: k = 0, controls the 1/z piece
  opts = odeset('RelTol', reltol, 'AbsTol', 1e-10, 'MStateDependence', 'none', 'Mass', @Mk);
  [tout, uout] = ode15s(@deqk, z, uin, opts);
  nk2s7 = uout(:, 1);
  nk2s7 = 2 + z'.*nk2s7; % nk2d with 1/z^2 factored out % checked feb 27

  save('nf.mat', 'nh2s7', 'nk2s7', 'na2t2', 'na2t4', 'na2t5');

  figure
  %plot(za, na2a, zb, na2b, zc, na2c, z, nh2b, z, nk2d);
  plot(z, na2t2, z, na2t4, z, na2t5, z, nh2s7, z, nk2s7);
  %plot(z, nh2b, z, nk2dm, z, nk2c);

end

